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Motivated by the experimental study of Tayebi et al. [Nature Mater. 11, 1074 (2012)] on phase 
separation of stacked multi-component lipid bilayers, we propose a model composed of stacked 
two-dimensional Ising spins. We study both its static and dynamical features using Monte Carlo 
simulations with Kawasaki spin exchange dynamics that conserves the order parameter. We show 
that at thermodynamical equilibrium, due to strong inter-layer correlations, the system forms a 
continuous columnar structure for any finite interaction across adjacent layers. Furthermore, the 
phase separation shows a faster dynamics as the inter-layer interaction is increased. This temporal 
behavior is mainly due to an effective deeper temperature quench because of the larger value of the 
critical temperature, T c , for larger inter-layer interaction. When the temperature ratio, T/T c , is 
kept fixed, the temporal growth exponent does not increase and even slightly decreases as function 
of the increased inter-layer interaction. 


I. INTRODUCTION 

Biological membranes are constructed out of two 
monolayers (leaflets) arranged in a back-to-back configu¬ 
ration. They are mainly composed of phospholipids but 
contain also other molecules such as cholesterol, glyco- 
sugars, and proteins pQ. In living organisms, these mem¬ 
branes can form multi-lamellar stacks known as lamel¬ 
lar bodies |2]. Examples of such highly folded membra¬ 
nous structures are thylakoid membranes of photosyn¬ 
thetic cyanobacteria or plant chloroplasts, and stratum 
corneum of human skin. Since multilamellar structures 
can combine single membrane functions in series, they 
offer possibilities for novel applications in photonics and 
as bio-sensors. 

Over the last decade, many studies have been per¬ 
formed on artificial giant unilamellar vesicles (GUVs) 
composed of ternary mixtures of saturated lipid such 
as sphingomyelin, unsaturated lipid such as DOPC (1,2- 
dioleoyl-sn-glycero-3-phosphocholine) and cholesterol [3, 
3]. By decreasing temperature, these ternary mix¬ 
tures undergo a lateral phase separation, where a liquid- 
disordered (Ld) phase coexists with a liquid-ordered (L 0 ) 
one. It is known that the L 0 phase is rich in saturated 
lipid and cholesterol, while the Ld phase is rich in the 
unsaturated lipid. 

In a recent experimental study, Tayebi et al. [5] re¬ 
ported that a stack (typically composed of several hun¬ 
dred layers) of multicomponent lipid bilayers with phase- 
separated domains exhibits inter-layer columnar order¬ 
ing. Using ternary mixtures of sphingomyelin, DOPC 
and cholesterol, it was observed that domains in stacked 
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bilayers align one on top of the other, thereby forming an 
uninterrupted columnar ordering across hundreds of bi¬ 
layer membranes. Such a cooperative multilayer epitaxy 
was attributed to the interplay between intra-layer do¬ 
main growth and inter-layer coupling. The formation of 
columnar structures in stacked bilayers is important be¬ 
cause it allows for electrical currents and transport pro¬ 
cesses to pass through many transmembrane channels in 
a cooperative and efficient manner. Other possible appli¬ 
cations of the columnar ordering can be as templates for 
membrane protein crystallization, which is necessary for 
X-ray structural analysis of membrane proteins incorpo¬ 
rated in bilayers. 

As far as the dynamics of phase separation in stacks 
of membranes is concerned, the temporal evolution of 
the average inplane domain size, R , was shown to obey 
a power-law growth, R ~ t a with a ~ 0.455 [5]. This 
exponent is larger than the value obtained using GUVs 
with a single bilayer, for which the reported experimental 
value is a « 0.28±0.05 [61. Hence, Tayebi et al. concluded 
that the inplane domain growth in each of the bilayers of 
the stack is faster, as compared to the domain growth in 
GUVs. 

In a subsequent paper [7] , a model based on regular so¬ 
lution theory, which takes into account the inter-lamellar 
coupling of inplane phase-separated domains, was pro¬ 
posed. The calculated phase diagram was presented in 
terms of intra-layer and inter-layer coupling parameters, 
and contains three different regions: (i) a “one-phase” 
region in which the system does not exhibit phase sepa¬ 
ration; (ii) a “two-phase” region in which two phases co¬ 
exist and domains in different layers along the normal z- 
direction are completely aligned and have the same com¬ 
position in the various layers, and (iii) a “multi-phase” 
region in which there are unaligned inplane domains with 
different composition in the different layers. According to 
Ref. 2], the transition line between the “two-phase” and 



2 


“multi-phase” regions strongly depends on the number 
of layers in the stack which was varied up to ten layers. 
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FIG. 1. (a) Schematic illustration of a stack of binary mem¬ 

branes, taken here as a stack of three bilayers in the 2 -direction. 
Each bilayer is composed of two identical leaflets containing sat¬ 
urated lipids (A, black) and unsaturated lipids (B, white). Sat¬ 
urated and unsaturated lipids typically form L 0 and La phases, 
respectively. As the lipid molecules are not allowed to exchange 
between different bilayers, their composition in each bilayer is 
fixed, (b) The stacked two-dimensional (2d) Ising model. Here 
the bilayer structure of each membrane is neglected. Lipid A 
and B correspond to spin up (black) and spin down (white), re¬ 
spectively. J is the coupling parameter between nearest-neighbor 
spins in the same layer, while J' is the coupling parameter be¬ 
tween spins belonging to two nearest-neighboring layers. 

Being motivated by these works 00, we investigate 
the correlation between lateral phase separation in a 
stack of multi-layer membranes using a spin model called 
the stacked two-dimensional (2d) Ising model. This is the 
simplest model to describe a stack of binary membranes 
composed of two types of lipids. The model is the same 
as the anisotropic three-dimensional (3d) Ising model for 
a finite stack in the 2 -direction. The important differ¬ 
ence between the two models is that in the former the 
order parameter (magnetization) in each layer is con¬ 
served. This requirement is based on the experimental 
fact that the A/B lipid composition in each layer almost 
does not change during experimental times. 

In our model, we study the thermodynamical equi¬ 
librium features using Monte Carlo (MC) simulations. 
The main reason that we performed MC simulations 


rather than analyzing the mean-held free energy describ¬ 
ing phase separation (as studied in Ref. 0 ), is to allow 
us to investigate the role of thermal fluctuations on the 
inter-layer domain correlation in stacked membranes. We 
show that the domains in each layer are correlated along 
the vertical 2 -direction, for any finite value of the inter¬ 
layer interaction is positive, i.e., J' > 0. Hence, the 
system is either in a one- or two-phase state in equilib¬ 
rium, and in our model the “multi-phase” state is not 
obtained in the thermodynamic limit of infinite lateral 
size, as long as the inter-layer coupling J' > 0. As antici¬ 
pated, it is found that the phase-transition temperature, 
T C ( J'), increases as function of the inter-layer interaction 
parameter. 

We also investigate the dynamics of phase separation 
at fixed temperature T in the two-phase coexistence re¬ 
gion. We show that the accelerated temporal behavior of 
the phase separation for the stack is mainly driven by the 
increase of the temperature quench, AT = T C {J') — T, 
because T c (,7') becomes larger for larger J'. However, 
if the ratio T/T C (J') is kept fixed, the dynamics of the 
phase separation is actually slower for larger values of the 
inter-layer coupling, J'. 

In the next section, we describe the stacked 2d Ising 


model and review the MC simulation method. In Sec. Ill 


we present the equilibrium properties of the model, and 
discuss the condition for domain columnar ordering. Sec¬ 
tion IV describes the dynamics of domain growth for dif¬ 
ferent values of the inter-layer interaction, and it is com¬ 
pared with a previous theoretical work. 


II. MODEL AND SIMULATION TECHNIQUE 


In our simulations, we use the stacked 2d Ising model , 
shown in Fig. [jja). We consider a stack of two- 
component lipid bilayer membranes composed of an A/B 
lipid mixture, although the experimental systems often 
consist of ternary lipid/cholesterol mixtures. This sim¬ 
plification does not affect the essential feature of the lat¬ 
eral phase separation. Another simplification is that we 
treat only symmetric bilayers where the composition of 
the two leaflets is identical. Hence, each lipid bilayer 
having a finite thickness can be mapped into a 2d Ising 
model with conserved magnetization, expressing the fact 
that no lipid is allowed to exchange across layers. The 2d 
Ising layers are stacked in the 2 -direction, and they inter¬ 
act with their two nearest-neighboring layers, as depicted 
in Fig.JTJb). 

The Hamiltonian of this stacked and coupled 2d Ising 
system can be written as: 


H = - J ^2 Si,pS itP ' - J' ^2 Si,pSi+ i, P 

i,(P,P') i >P 

— ^ ( 1 ) 


where up/down values of the spin, Si tP = ±1, at p = 
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FIG. 2. Time evolution of phase separated domains in the stacked 2d Ising model at different MC steps for A = 0.1 and T/J= 1.63. 
The other parameters are Si, p = 0 and L = L z = 64. For presentation purposes, only the interfaces of domain boundaries are 
shown, and the two different sides of the interfaces are represented by green and brown. The system is fully equilibrated after about 
10 7 MCS. 


(x, y) in the i -th layer corresponds to a lattice site oc¬ 
cupied by an A or B lipids, respectively. The coupling 
between nearest-neighbor spins in the ary-plane (denoted 
by ( p,p')) is J, while the coupling with the nearest- 
neighbor spins across layers in the z-direction is J'. The 
physical origin of the inter-layer interaction J' is primar¬ 
ily attributed to direct van der Waals attractive inter¬ 
actions acting between neighboring bilayers [5]. Other 
non-specific interactions, such as electrostatic and/or hy¬ 
dration interactions, can be taken into account through 
the second virial coefficient and will affect the value of J' 
as well [9l HO]. Throughout this paper, we shall use the 
dimensionless ratio defined by A = J' / J as a measure of 
the inter-layer coupling strength. 

In the above Hamiltonian, y,i is the external field 
(chemical potential), which fixes the average magneti¬ 
zation (A/B composition) in the i-th layer. Although 
Hi can, in general, take different values for different lay¬ 
ers, we consider here the case where all of them are the 
same, y,i = y, fixing the same value of lipid composi¬ 
tion in all layers. This assumption holds also for the 
dynamical states since we do not allow the lipids to be 
exchanged across different layers. The average order pa¬ 


rameter (A/B composition) in the *-th layer is denoted 
by Si,p, and throughout this paper (except in Fig. [BJb)) 
we choose S',. p = 0, which corresponds to a symmetric 
1:1 A/B lipid mixture, i.e., at the critical composition. 
This is equivalent to setting the value of the chemical 
potential to zero, i.e., y = 0. 

The present model is related to the anisotropic 3d Ising 
model for a finite slab. The special case of A = 1 corre¬ 
sponds to the isotropic 3d Ising model, whereas for A = 0 
the stack is composed of non-interacting 2d Ising layers. 
One interesting issue related to the anisotropic model, 
0 < A < 1, is the crossover from 2d to 3d critical behav¬ 
ior m that will be explored below. We also note that 
the stacked 2d Ising model has been studied a great deal 
in connection with multilayer adsorption phenomena on 
attractive substrates mis], but not in the context of 
layers of binary mixtures with conserved magnetization 
(order parameter) as studied in this paper. 

We investigate both the statics and dynamics of a stack 
of membranes based on the Hamiltonian presented in 
Eq. 0. We employ MC simulations for classical Ising 
spins on a finite L x L x L z lattice. Periodic bound¬ 
ary conditions are used in all three directions. The spin 
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configurations are updated using Kawasaki exchange dy¬ 
namics jl4l in order to conserve the magnetization in each 
layer. This is based on the experimental fact that lipids 
almost do not exchange across different layers. Hence, 
their A/B inplane composition is fixed during experimen¬ 
tal times. 


The MC simulations presented here are performed in 
the following way. At each MC trial step, a site on the 3d 
lattice and one of its nearest neighbors in the same layer 
are chosen at random. If the two spins are alike, a new 
site is again chosen at random. This process is repeated 
until two unlike nearest neighbor spins are found. Then, 
the probability of exchanging the two spins is determined 
by the standard Metropolis algorithm na. If the energy 
difference due to the spin exchange becomes negative, i.e., 
A E < 0, we accept the exchange. Otherwise, we accept 
the exchange with a probability exp(—A E/T), where T 
is the temperature and the Boltzmann constant, /cb, was 
set to unity. 


In one Monte Carlo step (MCS), this procedure is re¬ 
peated L x L x L z times. The MC simulations are car¬ 
ried out by annealing the temperature gradually from an 
initial infinite temperature for which the spin configura¬ 
tions are completely disordered and uncorrelated. The 
first 10 5 (or in some cases up to 10 6 ) MCS are discarded 
in order to reach thermal equilibration. Furthermore, to 
avoid correlations between different equilibrated configu¬ 
rations, measurements are taken every 20 MCS, and we 
averaged over 10 5 independent system configurations, in 
order to obtain sufficient statistics. 




In order to investigate the phase separation dynam¬ 
ics, we monitor the domain coarsening as a function of 
time (MCS) at a constant temperature below T c . An 
example of a typical time evolution of phase separation 
is presented in Fig. [2] for A = 0.1, T/J = 1.63 and 
L = L z = 64, where six snap-shots are shown from 
10 2 MCS till 10 7 MCS. For clarity purposes, only the 
boundaries between domains of spin up (rich in lipid A) 
and spin down (rich in lipid B) are shown. In the ini¬ 
tial time steps, the phase separation occurs inplane, and 
the domains coarsen without much out-of-plane coupling 
(due to the rather small value of A = 0.1). As time 
evolves, the inplane coarsening is also followed by out-of- 
plane columnar ordering, where the lipid A (and lipid B) 
rich domains are highly correlated along the z-direction. 
This is clearly seen for the fully equilibrated configuration 
occurring after about 10 7 MCS (last snap-shot). Here the 
two color boundaries, represent the two sides of the do¬ 
main boundaries (while the inside of the domain is not 
shown). The boundaries look like extended interfaces 
separating inplane domains that are vertically connected 
along the z-direction, in agreement with experiment [5]. 


FIG. 3. (a) Specific heat per lattice site, c, as function of 

the dimensionless temperature T/J, for different lateral system- 
size L = 16,24,32,40,48. The other parameters are Si, p = 0, 
A = 0.1 and L z = 8. For each system size, the peak position of c 
is identified with an effective “phase transition" temperature, (b) 
Finite-size scaling analysis of the phase-transition temperature, 
T C (L)/J for A = 0.1. The apparent phase-transition tempera¬ 
ture is plotted as a function of 1/L. The solid line is the fit given 
by Eq. ([3]) with v = 1 (see text). The extrapolated value for the 
critical temperature is T c (A = 0.1)/J = 2.85. 


III. STATIC PROPERTIES OF THE STACKED 
DOMAINS 


In order to determine the phase-transition temperature 
and obtain the corresponding phase diagram, we compute 
the specific heat per lattice site defined as 


c = 


1 1 
L 2 L Z T 2 


(( H 2 > - (IT) 2 ) , 


( 2 ) 
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where H is given by Eq. ([I]) and (■ ■ •) indicates an en¬ 
semble average. We note again that the above specific 
heat is calculated at constant magnetization (correspond¬ 
ing to constant lipid concentration in our model) of each 
layer. In our simulations, the ensemble average is taken 
by averaging over independent equilibrium spin configu¬ 
rations as explained in Sec. [TTJ For a given system size 
and dimensionless ratio A, we calculate c as function of 
the dimensionless temperature T/J. Such a dependence 
of c on T/J is presented in Fig. [3ja) for several lateral 
system-sizes, L, and for A = 0.1, L z = 8, recalling that 
L z is the number of layers of the 3d stack. 

For each system size, we associate the peak position of 
the specific heat with the apparent critical temperature, 
T C (L, A), for a system with a finite size, L. Finite-size 
scaling analysis is then performed in order to determine 
the critical temperature for a slab of a finite L z layers 
in the thermodynamic limit (L —> oo). In Fig. [3](b), we 
plot T C (L, A = 0.1) as a function of 1/L for the same 
parameters as in (a). The plotted data are fitted with 
the following finite-size scaling assumption: 

T c (L,X) = T c (X) + aL~ 1 / 1 ', (3) 

where T c (A) = T C (L —> oo, A) is the infinite system criti¬ 
cal temperature for a given A, a is a non-universal prefac¬ 
tor, and v is the 2d critical exponent for the correlation 
length in the xy-plane. We set v = 1 in our analysis, fol¬ 
lowing the work by Pham Phu et al. m , who performed 
extensive MC simulations on magnetic Ising films (with 
A = 1) PU- We choose this 2d critical exponent for the 
fitting because it was shown m that the 2d character of 
the film is dominant even for L z = 13. The extrapolated 
critical temperature for L —> oo obtained from Fig. [3](b) 
is T c (X = 0.1 )/J = 2.85. We repeat this procedure for 
different values of the inter-layer interaction parameter in 
the range of 0 < A < 1, and determine the corresponding 
critical temperature, T c (A). We note that the value v = 1 
provides a good fitting for all the A values examined. 

Somewhat surprisingly, finite-size effects in the z- 
direction are much weaker as compared to those in the 
lateral direction. This is shown in Fig. [4j where we plot 
c as a function of T/J when (a) A = 0.1, (b) A = 0.5 and 
(c) A = 1 for different number of layers, L z = 4,8,12,16, 
while the lateral size L = 48 is kept fixed. For all A 
values studied here (0.1 < A < 1), the observed peak 
position: T/J ss 2.65 in (a), 3.30 in (b) and 4.10 in (c), 
is almost independent of L z , at least for L z > 8. This 
means that, in our model with a fixed imposed magneti¬ 
zation (A/B composition) in each layer, the correlation 
in the ^-direction is very strong due to the cooperative 
behavior of domains in different layers. 

For fully equilibrated configurations, as shown in Fig. [2] 
after 10 7 MCS, the domains are highly connected verti¬ 
cally along the ^-direction, from the bottom layer to the 
top one. This is also shown in Fig. [5] in which the colum¬ 
nar structure of domains in different layers is clearly 
shown. Hence, the correlation length in this direction ex¬ 
ceeds L z , and the constraint of fixed magnetization (A/B 


composition) in each layer induces a strong structural 
correlation in the ^-direction even though the inter-layer 
interaction J' is smaller than the intra-layer interaction 
J (A < 1). A more quantitative argument for the domain 
connectivity will be given later. Because the number 
of layers, L z , barely affects the MC results as shown in 
Fig.[4]for A = 0.1, 0.5 and 1, most of the simulations were 
done using L z = 8, which is sufficiently large in our case 
to observe the asymptotic behavior of L z —> oo. For the 
anisotropic 3d Ising model without any constraint of con¬ 
served magnetization, as previously studied in Ref. m, 
a very weak system-size dependence of the apparent crit¬ 
ical temperature was observed by measuring the planar 
susceptibility. 

The results of finite-size scaling analysis are shown 
in Fig. [6j where we plot T c as a function of A. The 
critical temperature interpolates between the 2d and 3d 
Ising results, T 2d < T C (A) < T ( ? d ; the exact value in 
2d (corresponding to A = 0) is known to be T 2d /J = 
2/ln(l + a/2) ~ 2.269 for square lattices fTS], and the 
numerical estimate in 3d (corresponding to A = 1) is 
T dd / J ss 4.511 for cubic lattices [T5] , These two limits 
are recovered in our simulations and are seen in Fig. [6] 
for A = 0 and 1, respectively. Although a more detailed 
A-dependent scaling behavior of T C (A) was previously dis¬ 
cussed in the limit of very small A puns, we shall gen¬ 
eralize the argument for the anisotropic case of finite A, 
0 < A < 1. When T < T c (A), the stack undergoes a 
phase separation, and the inplane domains rich in lipid A 
(spin up) are interconnected along the ^-direction, bridg¬ 
ing between adjacent layers and forming large connected 
domains of the same average composition. The same fea¬ 
ture also occurs for the B-rich domains. Such a behavior 
can be clearly observed in Fig. [5} 

In order to monitor quantitatively the degree of inter¬ 
connectivity of domains in different layers, we define the 
following quantity: 

« 2 = /(i;(/X>.p-3./)). (4) 

where the average is taken over equilibrated MC config¬ 
urations as explained above. This quantity can be cast 
also as: 

^ = L 2 L 2 ~ S iiP )(S jtP — Sj tP )) , (5) 

z P i,j 

and represents a special “magnetic susceptibility”, where 
the correlations are taken only along the ^-direction and 
then averaged laterally in each of the planes. When the 
domains are connected along the z-direction, the summa¬ 
tion over different i-layers will produce a large value of 
6, while S is small if the domains are uncorrelated across 
the layers even for T < T c ( A). In Fig. [7J we plot 6 2 as 
a function of T/ J for different values of A, while fixing 
L — 16 and L z = 8. Notice that even for A as small as 
0.05 (blue diamonds), S 2 tends to increase as the tem¬ 
perature decreases below T c (A). This means that the 
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FIG. 4. Specific heat per lattice site, c, as function of the dimensionless temperature T / J , for different systems size L z = 4,8,12,16 
for (a) A = 0.1, (b) A = 0.5 and (c) A = 1. The other parameters are Si tP = 0 and L = 48. The observed peak position: T/J m 2.65 
in (a), 3.30 in (b) and 4.10 in (c), is almost independent of L z , at least for L z > 8. 



FIG. 5. Time evolution of phase-separated domains in a stacked 2d Ising model of eight layers, L z = 8, at different MC steps for 
(a) Si tP = 0 and (b) Si yP = 0.4. The other parameters are A = 0.1, T/J = 2.0 and L = 256. 


domains are connected in the ^-direction once the phase 
separation takes place. On the other hand, domains are 
independent and uncorrelated only when the inter-layer 
interaction is extremely small, i.e, A < 0.001 in Fig. [7] 
The situation is found to be marginal when A = 0.01 (red 
triangles) because S 2 then slightly deviates from zero at 
low temperatures. 

Based on our MC results, we conclude that in the ther¬ 
modynamic limit, L —> oo, domains will always be con¬ 
nected for any finite inter-layer interaction, J' > 0. We 
give now a simple argument supporting this conclusion, 
and show that in the limit L —> oo but with a finite 
number of layers, L Zl the domains in different layers are 
uncorrelated only when J' = 0 (A = 0) is strictly obeyed. 
For the symmetric A/B case (SYp = 0), each layer will 
eventually phase separate into two semi-infinite domains: 
one composed by the A lipid (spin up) and the other by 
the B lipid (spin down), as shown schematically in Fig. [8] 
When the domains are fully correlated in the z-direction, 


as in Fig. [8](a), the total free energy of the stack consists 
of the contributions: 

-Peon = — J' L Z L 2 + -Pjntra, (6) 

where Fi ntra accounts for the intra-layer interactions. On 
the other hand, when the inplane domains are completely 
random and disconnected, as sketched in Fig. |8](b), the 
total free energy is dominated by an entropy contribution 
of arranging a random stack of A and B domains along 
the z-direction, 


Pdis — —TL Z In 2 + idntra, (7) 


with the same Flntra as before because this term is com¬ 
mon for both free energies. By comparing Eqs. © and 
([T]), the threshold inter-layer interaction, (J 1 )*, separat¬ 
ing the two states, is given by: 


(j'r 


T In 2 
L 2 


( 8 ) 
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FIG. 6. The phase-separation temperature, T c (\)/J, at the 
critical composition, as a function of the interaction parameter 
A for symmetric A/B mixtures, Si, p = 0. The system is in a 
phase-separated state below the solid line, and in a one-phase 
state above the line. 


Notice that ( J')* depends on L but not on L z . For finite 
temperatures, it vanishes in the thermodynamic limit of 
L —f oo. Hence, this simple scaling argument suggests 
that domains are always connected in the 2 -direction for 
any finite value of J'. Therefore, for all A > 0, in the 
phase-separated region (below the critical temperature) 
presented in Fig. [6j domains should always form inter¬ 
connected structures along the 2 -direction. As shown in 
Eqs. ([6]) and (J7|, the internal energy scales with L 2 , while 
the entropy due to the random stacking of domains does 
not depend on L. Hence, the entropic effect can never 
overcome the internal energy in the thermodynamic limit, 
and leads to the stability of the columnar structure. This 
conclusion is not in agreement to that of Tayebi et al. if], 
who claimed that there is a “multi-phase” state in which 
domains are not aligned and have different compositions 
even in thermodynamical equilibrium. 

In the simulations, (J')* can be finite due to finite-size 
effects. For instance, if the temperature is chosen to be 
T/J = 1 in Fig. [?J the threshold value for L = 16 can 
be estimated as A* = « 2.7 x 10 -3 . Since A = 

10~ 2 (red triangles in Fig. f7[) exceeds this threshold, the 
corresponding S 2 takes larger values at low temperatures. 
Moreover, the very weak finite-size effects along the 21 - 
direction is consistent with the lack of L z -dependence of 
(J')* in Eq. 


IV. DYNAMICS OF PHASE SEPARATION 

We address now the effects of inter-layer interaction 
on the dynamics of phase separation as the system con- 



FIG. 7. The out-of-plane domain connectivity, S 2 , defined in 
Eq. 0, as a function of the dimensionless temperature T/J, for 
different values of A = 0, 0.001, 0.01, 0.05, 0.1, 0.5, 1.0. The 
other parameters are Si, p = 0, L = 16 and L z = 8. The transi¬ 
tion temperatures for different A values are indicated by arrows. 
The value of S 2 becomes larger when domains are correlated 
along the 2 -direction between different layers. This increase in 
S 2 is observed for lower temperatures and larger A. 



FIG. 8. Schematic representation of phase separated domains 
in a stack of membranes. Black and white domains are rich in 
A and B lipids, respectively. Two extreme cases are shown; (a) 
domains are fully connected in the 2 -direction, (b) domains are 
arranged at random and are disconnected. 

verges towards its thermal equilibrium state. Under the 
assumption that scaling laws can be applied, the average 
domain size R increases according to a temporal power- 
law: R(t) ~ t a 0. For 2d systems for which the total 
domain area is conserved, the average domain size R is 
inversely proportional to the total interface length £, i.e., 
R ~ £ _1 [21] [22'. This can easily be seen because R 
and £ are related by £ = 2irnR ~ nR , where n is the 
number of domains, and the total area of all domains, 
A = ttuR’ 2 ~ ni? 2 , is a conserved quantity. Hence, within 
the scaling hypothesis, the total interface length (in 2d) 
should behave as 

£{t) ~ t~ a . (9) 

In our stacked Ising model, we calculate the interface 
length in each of the layers and average it over different 
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FIG. 9. (a) The temporal evolution of the total interface length 
l as a function of time (MCS) for different values of A = 0, 0.2, 
0.4, 0.6, 0.8, 1.0, and for a temperature quench from the one- 
phase state (T — \ oo) into the two-phase state at T/J = 2.0. 
The A/B mixture is symmetric, Si tP = 0, L = 256 and L z = 8. 
The average over three independent MC runs is taken for each 
A value. The two dashed lines represent a power-law behavior 
with exponent a = 0.07 and 0.28, which roughly bound the 
two limiting behaviors of the A-dependent exponent, a. (b) The 
domain growth exponent a as a function of A, as obtained from 

(a)- 


layers. Note that the total interface length is propor¬ 
tional to the first term of the Hamiltonian in Eq. ([!]), 
which enumerates the number of bonds across the in¬ 
plane A/B interface. 

In Fig. [9^a), we plot the temporal evolution of the 
total interface length in 2d, l(t), (and averaged along 
the ^-direction), as a function of time measured in MC 
steps. The temperature quench into the two-phase region 
is done for a fixed temperature, T/J = 2.0 < T c (A), in 



MCS 



X 


FIG. 10. The temporal evolution of the total interface length 
l as a function of time (MCS) for different values of A = 0, 
0.2, 0.4, 0.6, 0.8, 1.0, and for a temperature quench from the 
one-phase state into the two-phase one, with final temperature 
satisfying T = 0.6T C (A). The A/B mixture is symmetric, Si, p = 
0, L = 256 and L z = 8. The average over three independent MC 
runs is taken for each A value. The two dashed lines represent a 
power-law behavior with exponents a = 0.14 and a = 0.24. (b) 
The domain growth exponent cr as a function of A, as obtained 
from (a). 


order to mimic the experiment that is conducted at fixed 
room temperature. Several values of A are studied, and 
the other parameters are L = 256 and L z = 8, with av¬ 
erages taken over three independent MC runs. For each 
A value, the scaling behavior of Eq. ([9]) is analyzed, and 
we extract the growth exponent a from the late stage 
kinetics. We find that for A = 0 (2d case), the growth 
exponent has the smallest value of a « 0.07, while for 
A > 0, it is a function of A and increases up to a ss 0.28, 
as shown in Fig. |9](b). 
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Although this result may explain the fact that the 
phase separation has an accelerated dynamics in stacked 
membranes as compared to GUVs (isolated single mem¬ 
branes), we should keep in mind that T C (A) increases as 
function of the inter-layer coupling A > 0, as shown in 
Fig-! As long as the final quench temperature is fixed 
to T/J = 2.0, the temperature quench depth defined by 
AT = T c (A) — T becomes larger as the value of A is in¬ 
creased. This may explain why the growth exponent a 
becomes larger with increasing A, for a fixed T-quench. 

In order to have a better comparison between differ¬ 
ent A values, we evaluate in Fig. [To] the growth exponent 
in a different way. We now keep a constant quench ra¬ 
tio T/T c { A) = 0.6, where T is the final quench tempera¬ 
ture, and the critical temperature T c (A) depends on A, as 
shown in Fig. [6j For these deeper temperature quenches 
(farther from T c (A)), the estimated growth exponent is 
a « 0.24 for A = 0 (pure 2d case), and 0.13 < a < 0.16 
for 0.2 < A < 1.0. Note that the a-values are only weakly 
dependent on A > 0. 

Finally, we elaborate on the decreasing A-dependence 
of the growth exponent a, and show that this behav¬ 
ior is consistent with the change in the dimensionality 
of the stack from 2d to 3d. In general, the growth ex¬ 
ponent associated with phase separation depends on the 
dimensionality | 23 | . In this context, we mention the scal¬ 
ing argument of Binder and Stauffer on phase-separation 
dynamics of particles that undergo cluster reaction and 
diffusion processes [53j. Under the assumption that most 
particles that leave a cluster reimpinge on the same clus¬ 
ter at later times, the diffusion coefficient D of a cluster 
of size R was shown to scale as D ~ R~^ 1+d \ where d 
is the embedded space dimension. If we further assume 
that the domain size R is the only length scale in the 
system, the scaling relation for a simple diffusion process 
is given by R 2 ~ Dt. This argument yields the growth 
exponent to be a = 1/(3 + d). Hence, the predicted val¬ 
ues from this scaling conjecture are a = 1/5 for d = 2 
and a = 1/6 for d = 3. 

Our simulation results, namely, a « 0.24 for A = 0 
and a « 0.14 for A > 0.2 compare favorably with this 
prediction. The growth exponent decreases for finite A 
because the system crosses-over from 2d to 3d. This is 
due to the fact that the growing phase-separated domains 
are inter-connected along the ^-direction for A > 0. It 
should be noted, however, that the absolute value of a 
obtained from the simulation is not universal but strongly 
depends on the quench depth as shown in Fig. [9] This 
explains why the above exponents are not in complete 
agreement with the simple scaling argument of Binder 
and Stauffer. 


V. CONCLUDING REMARKS 

Motivated by recent works of Tayebi et al. 0E3, who 
studied experimentally and theoretically the phase sep¬ 
aration in stacks of multi-component lipid bilayers, we 


have investigated the stacked 2d Ising model given in 
Eq. 0. We use a Monte Carlo simulation scheme with 
Kawasaki exchange dynamics that conserves the order 
parameter in each layer, in order to investigate both 
equilibrium and dynamical features. Performing finite- 
size scaling analysis only in the lateral direction, while 
keeping the stack thickness fixed (mimicking the exper¬ 
iment), we determine the phase-transition temperature, 
T C (A), by changing the inter-layer interaction parameter 
A = J'/ J. As shown in Fig. |6j the phase-transition tem¬ 
perature interpolates between that of the 2d and 3d Ising 
model. 

One of our main conclusions is that domains in each 
one of the layers are always interconnected along the z- 
direction, forming a continuous columnar structure for 
any finite inter-layer interaction J' > 0, as shown in 
Fig-11 This domain structure is in accord with the exper¬ 
imental findings for stacks of few dozen to few hundred 
layers [5]. However, the “multi-phase” region in which 
there are unaligned inplane domains with different com¬ 
position, as was predicted in Ref. [7j, is not found in 
our study at thermal equilibrium. Of course that such 
a “multi-phase” state can be transiently observed before 
the system reaches its fully equilibrated state, as can be 
observed in Figs. [2] and [5] 

We have also investigated the temporal evolution of 
domain formation in the stacked 2d Ising model. When 
the inter-layer interaction A increases, the phase sepa¬ 
ration appears to have an accelerated dynamics as can 
be seen by the larger values of the growth exponent, a, 
shown in Fig. 0b). However, these larger a values are 
mainly due to an increase in the phase-transition tem¬ 
perature, T c (A), as function of A; thus, a larger effec¬ 
tive temperature quench, AT = T C (A) — T, for fixed T. 
When the final temperature quench T is fixed relative 
to the phase-transition temperature as shown in Fig. [0] 
for T = 0.6T C (A), the growth exponent even decreases 
as the A value is increased. Our numerical findings for 
the growth exponent a are different than the value of 
a « 0.455, as found in the experiment [5]. One possi¬ 
ble explanation for this discrepancy can be the lack of 
hydrodynamic interactions in our MC simulations [5]. 

In this work, we have mainly discussed the case of 
Si, P = 0, corresponding to the critical composition of 
the A/B lipid mixture. Currently, we are investigating 
the dynamics of phase separation for off-critical compo¬ 
sitions, Si^p ^ 0 [see Fig. 0b)]. For such compositions, 
the phase-transition temperature is smaller than the crit¬ 
ical temperature. In the present simulations, the aver¬ 
age A/B lipid composition (order parameter of the Ising 
model) in each bilayer is restricted to stay the same. In 
the future, we plan to study membrane stacks where each 
layer has a different but fixed composition |2j3j. Further¬ 
more, since it is known from simulations that the pres¬ 
ence of a supporting solid substrate affects the dynamics 
of membrane domain growth |25J, it will be of interest to 
incorporate this substrate effect in future studies. 
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